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ABSTRACT 

Purely hydrodynamic numerical experiments into the evolution of astrophysical discs typi- 
cally include some sort of viscosity in order to cause accretion. In this paper, we demonstrate 
an alternative method of implementing viscous forces, with extremely good angular momen- 
tum conservation properties. The method is based on altering the cell fluxes, rather than incor- 
porating a viscous force. We test this method on the classical 'ring spreading' problem, and 
demonstrate angular momentum conservation at the 10 -8 level. 
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1 INTRODUCTION 

Gas dynamics dominates the physics of many areas of astronomy. 
Unfortunately, gas flow is a complicated subject, which has proven 
stubbornly resilient to analytic study. This is due to the non-linear 
nature of the Navier-Stokes equations which describe the flow. We 
are forced to look to computers to conduct numerical experiments, 
in order to better understand the cosmos. The reason that computa- 
tional work is better described as a 'numerical experiment' rather 
than a 'simulation' is the plethora of techniques available, with no 
obviously 'right' met hod. A recent compa rison of different codes 
on is the work of lde Val-Borro et"al]j200tj|) . who conducted a com- 
parison of a variety of codes on the planet-disc interaction problem. 
In this paper, we shall introduce an alternative method of imple- 
menting viscosity in a hydrodynamics code, and test it in the case 
of an accretion disc. 

Viscous processes in astrophysics are typically thought to 
originate from magnetic ef fects (e.g. the Magnetorotational Insta- 
bility (MRI) of Balbus and Ha wldll99ll is thought to dominate 
transport in most accretion discs), rather than purely hydrodynamic 
processes (which would be far too weak, given the low densities in- 
volved). However, magnetohydrodynamic (MHD) calculations are 
computationally expensive. It is therefore common to incorporate 
conventional viscous terms into the equations used, to approximate 
MHD effects. Observational effects of accretion can thereby be 
studied, without the added complexity of a full MHD calculation. 
However, since the details of MHD transport will certainly be dif- 
ferent from a simple physical viscosity, MHD calculations are re- 
quired to verify results fr om the physical viscosity approach. This 
has been demonstrated by Nelson (2005) in the context of the mi- 
gration of low mass planets, and by Wi nters et alj ( 2003 ) for planets 
which can open a gap in a circumstellar disc. A further reason to 
incorporate viscosity into a hydrodynamics code is to ensure that 
results are not dominated by details of the algorithm. All codes 



for solving the equations of hydrodynamics involve some numeri- 
cal dissipation. This dissipation is dependent on both the algorithm 
used, and the resolution computed, and there is no general method 
of calculating what effect it has. For this reason it is best not to 
refer to the dissipation as 'numerical viscosity,' since the effect is 
unlikely to be precisely that of a diffusion equation. By introducing 
some physical viscosity into the hydrodynamic equations, we can 
hope to overpower the unknown dissipation a parameter we can 
control. 1 Of course, if too much physical viscosity is required to 
achieve this, then the quality of our numerical method is called into 
question. 

The structure of this paper is as follows: in section [2] we in- 
troduce the equations of hydrodynamics, and discuss how viscosity 
has been implemented previously in section|3| Our new approach to 
incorporating viscous terms is summarised in section|4| We discuss 
implementation details in section|5|and our tests in section|6| 



2 THE HYDRODYNAMIC EQUATIONS 

The equations of hydrodynamics have been known for several cen- 
turies, and may be written in a variety of forms. A common one is 
the following: 

g+V-(pv) = (1) 

5 + (v-V)v = -ivp-V* + V-Ty (2) 
at p 

where p is the density, v the velocity of the fluid, p the pressure, $ 
the gravitational potential and Tg is the viscou s stress tensor (see, 
e.g. lLandau and Lifshitzll 9591 : iBatchelorl 1967ft . The first equation 
describes the conservation of mass and the second, conservation of 
momentum. An equation of state closes the system of equations, 
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1 References to 'numerical viscosity' are better read as 'minimum level of 
physical viscosity required to dominate intrinsic numerical dissipation' 
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and additional terms may be added as required. The viscous stress 
tensor is written as 



v Vv+(Vv) T + U--=v (V-v)I (3) 



where v is the kinematic shear viscosity, and ( is the kinematic bulk 
viscosity. 2 The first term of Equation[3]represents the resistance of 
the fluid to shear forces, while the second measures the resistance 
to dilation. Note that the first term is symmetric, which ensures that 
solid body rotation does not give rise to viscous forces. 



3 PREVIOUS WORK 

Computer codes designed to study astrophysical fluids have in- 
cluded viscosity for a number of years. We shall review the meth- 
ods used in this section, starting with a brief review of how such 
codes work. This discussion concentrates on time-explicit grid- 
based codes. Other types of code exist, notably Smo othed Parti- 
cle Hydrodynamics (see, e.g. iBenzll 199(1 Monagharj |T99l for a 
description), but the new method for computing viscosity we are 
introducing here is not appropriate for them. The following dis- 
cussion is far from comprehensive (to be so would require several 
books), but serves as a brief outline of the principles involved. 

As their name implies, grid-based codes impose a grid on the 
computational volume, breaking it down into cells. The various 
flow variables (density, velocity, etc.) are stored for each cell. Sim- 
ulation time is advanced in a two step process. In the source or 
flux step, the fluxes of mass, momentum etc. through each cell face 
are computed. These fluxes are then used to update the cell quanti- 
ties during the transport step. This approach enforces conservation, 
since the flux out of one cell will be the flux in to another. 

Codes are further classified according to how they obtain 
their fluxes. Some do this by direct differencing of the equa- 
tions of hydrodynamics, and are often refe rred to as "ZE US-like" 
- a reference to the ZEUS code of IStone and Normarl |l992). 3 
The alternative method in common use is due to Godunov, which 
solves a 1-D shock tube problem at every cell interface (the so- 
called Riemann problem). Godunov's scheme is most commonly 
implement ed using th e Piec ewise Parabolic Method (PPM) of 
Colell aand Woodward 11984) . All explicit schemes for numerical 
hydrodynamics are subject to the Courant-Friedrichs-Lewy (CFL) 
condition for their maximum timestep. This requires that 

At < (4) 

\v x \ + c s 

Physically, this means that information must not propagate more 
than one grid cell per timestep. 

In codes like these, viscosity has generally been implemented 
in one of two ways. The first is to incorporate the viscous forces 
directly into the computation of the fluxes. The acceleration of the 
fluid is calculated as 



P 



Vp - V$ + V • Ti-j 



(5) 



where the viscosity has been added as an extra term to the usual 
forces of gravity and pressure gradient. An alternative approach 

2 If Equation[3]is written in terms of momentum, then these coefficients 
must be multiplied by the density, and r\ = pv is sometimes used 

3 Note that this does not imply that all these codes are derived from ZEUS; 
simply that ZEUS contains the most widely known implementation of the 
principles involved 



adds a separate substep to the source and transport steps outlined 
above. In this, the velocities are evolved separately, according to 



<9v _ m 



(6) 



This method is used in the FARGO code (written bv lMassetl2 00u). 
and also in ZEUS bv lStoneetai] fi999). In both cases, the viscous 
terms are being treated as a force. When used explicitly, both of 
these methods lead to an extra CFL condition of the form 



At < 



(Ax) 
2v 



(7) 



which must also be satisfied everywhere on the g rid. Since this can 
become somewhat restrictive, Kiev 1 1989, 1999) details an implicit 
solution method, which is not subject to Equation^ 



4 A NEW APPROACH 

In this section, we shall discuss a new approach to incorporating 
viscosity into hydrodynamics codes. 

Consider Equation|2| We can recast this in conservative form 

as 



^ + V • pvv = -Vp - V$ + V ■ Tij 
at 



(8) 



The pvv term is the momentum flux, which is supplied by the 
source step. Since the divergence is taken of both this and the 
viscous stress tensor, an alternative formulation is suggested: take 
the momentum fluxes, and subtract the viscous stress tensor. Equa- 
tion|S|then becomes 



^ + V ■ (pvv -T ij ) = -Vp - V$ 



(9) 



Instead of including viscosity into the source step, we should use 
it to generate an extra momentum flux, given by the viscous stress 
tensor itself. The transport step of the code will then act on the 
hydrodynamic and viscous fluxes. This has the advantage of being 
conservative, since whatever is added to one cell will be removed 
from another. We shall see how this leads to excellent conservation 
of angular momentum. However, since this method is still time- 
explicit, it remains subject to Equation^ As an additional bonus, 
this approach only requires numerical first derivatives, rather than 
the two derivatives required by the approach of Equations|5|and|6| 



5 IMPLEMENTATION 

In this section, we shall discuss the implementation of the approach 
outlined in section |4] in a real cod e. The code we shall use the 
FLASH code of lFrvxell et alJ l2OO0t) . an adaptive mesh refinement 
(AMR) code based around a PPM hydrodynamics solver. 4 

The downloadable source code contains a partial implementa- 
tion of the method outlined in section|4| (indeed, it was this which 
inspired the present work). However, this portion of the code is not 
fully implemented even in cartesians. When performing numerical 
experiments into the dynamics of accretion discs, polar geometry 
is highly desirable. This is because hydrodynamics codes generally 
conserve quantities along the grid axes. Hence, to conserve angu- 
lar momentum (rather desirable for an accretion disc), a polar grid 



The source code is available at |http : / / flash . uchicago ■ edu/| 



An alternative approach to viscosity 3 



is (almost) always required. The downloadable source doesn't con- 
serve angular momentum, and some extra modifications (not de- 
tailed here) were necessary to fix this problem. We take this angular 
momentum conserving code as the basis of our modifications. 

To compute the necessary fluxes, we must evaluate Equation|3| 
in polar co-ordinates. This an unpleasant exercise, but is built on 
standard results. We find that the six independent components of 



Ti-j are: 
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8v r 



(10) 



(11) 



(12) 
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2v 
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~dz~ J 
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(V-v) (13) 



(14) 



(V- 



where 



_ 1 drv 

V ■ v = - 



1 dv£ dv z 

dr r 86 dz 



(15) 



(16) 



From this point on, we shall concentrate on the 2D non- 
axisymmetric (r, 6) case. We therefore take z — v z = d z — 0, 
which results in a considerable simplification of these equations. 

At this point, numerical details intrude. In FLASH, all flow 
quantities (density, velocity etc.) are stored at cell centres. We wish 
to compute fluxes on the cell faces, and this introduces complica- 
tions. In particular, although T r< f, — T^ r mathematically, we wish 
to evaluate them in different places. We require T r< j, (that is, the 
flux of v r in the 6 direction) to be evaluated on the centre of the 6 
faces, while T^ r (the flux of in the r direction) is required on 
the centre of the r faces. This is illustrated in FigureQ 

We use centred derivatives, in order to achieve second order 
accuracy. For certain terms, this is easy. Consider the first term of 



t (f>r • 



dr 



-1,3 



(17) 



The second term is more troublesome. We want the derivative in 
the 6 direction, but centred on an r face. We solve this problem 
by evaluating the derivative at the cell centres on either side of the 
interface, and then taking the mean: 



1 dv r 
r 86 



h + Q 2 



n-i + r, 



where 



Q2 



fj+1 



„i,3+l 



~ 93-1 
„i,3-l 



9j+i ~9i-i 
The final term is comparatively straightforward: 



Vtf, 

r I 



i— 1,7 . i.i 
Vj. + Vj J 



"4> 



(18) 

(19) 
(20) 

(21) 



ri-x+n 

where the factors of two in the averages cancel top and bottom. 



Variables stored at cell centres 



IV 



i-1/2 



Figure 1. Storage of variables in FLASH. All variables are stored at cell 
centres, located at {(r^, rpj), (r^i, 6 A . . .}. The cell faces are located at 
the 'half positions. We require the components of the stress tensor to be 
centred on the cell faces 



Similar considerations apply when evaluating T r( p, although it is 
the first term which is trickier to centre correctly. Note that the dif- 
ferencing formulae given here break down close to r = (formally, 
once Ar w r). It would be possible to remedy this, but we shall 
not worry about the problem here. 

Since it is actually momentum density which is transported, 
all the computed fluxes must be multiplied by the density at the 
interface. This is computed as the mean of the density in the cells 
sharing the interface. Conversion to angular momentum is handled 
elsewhere in the code, and so does not affect the expressions given 
here. Boundary conditions are quite simple: no extra flux is added 
to faces on the edges of the computational volume. 



6 TESTS 

In this section, we shall discuss some tests of the new approach to 
implementing viscosity discussed above. 



6.1 The Viscously Spreading Ring 

The canonical test of viscosity in a computer c ode is arguably the 
viscously spreading ring described by Pringle 1 1981). In this test, 
the surface density evolution of material in a constant viscosity Ke- 
plerian accretion disc is monitored. If the surface density, E, is ini- 
tially a (5-function at r = ro, then at later times it is described by 



S(a;, t) oc j- ■ exp 



(22) 



where a; = r/ro,r = 12^rg and/j_ is a modified Bessel function. 

4 

The constant of proportionality is set by the mass of the ring. For 
further details, see lPringlel ll98ll) . Note that this only tests the T r ^> 
and T^ r components of the viscous stress tensor. Since 5-functions 
involve infinities, which are generally accepted to be numerically 
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Figure 2. Comparison of viscous ring evolution to analytic solution. Solid 
lines are results from FLASH, while dotted lines show the analytic solution. 
The r values corresponding to each pair of lines is also marked 



(a) Mass conservation as a function of r 



difficult, we start our tests from r > 0, and follow the subsequent 
evolution. 

To implement the spreading ring in FLASH, we use a gas with 
7 = 1.01 (since the Rieman solver cannot cope with an isother- 
mal equation of state). We reset the internal energy of the gas every 
timestep to maintain an aspect ratio h/r = 0.01, which is small 
enough to follow the analytic solution without introducing numeri- 
cal instabilities. We use a grid with n r — = 128 evenly spaced 
grid cells, with 0.5 < x = r/ra < 2 for the radial range. The 
inner and outer radial boundaries are reflecting, to enforce conser- 
vation. We initialise the surface density according to Equation 1221 
with r = 0.016. 

In Figure |2| we compare the evolution of the ring to the the- 
oretical solution. Although there are some deviations, particularly 
once the ring encounters the inner boundary (which was reflecting, 
not open), the agreement is excellent. Figure [5] demonstrates the 
conservation properties of the new method. Mass is conserved to 
the 10 -13 level (that is, machine precision), as one would expect. 
The rapid oscillations are obviously rounding error from one step 
to the next. Conservation of angular momentum is less good, but 
is still at the 10 -8 level. Computation of the viscous fluxes took 
approximately 2.5% of the total CPU time. In the absence of vis- 
cosity, angular momentum is conserved to the 10~ 13 level, and the 
ring does not spread. 

From this test, we see that the new scheme is performing ex- 
tremely well. 

6.2 Other Components 

The viscously spreading ring only tests the T r< j, component of the 
stress tensor. Unfortunately, we do not know of any generally ac- 
cepted tests of the other components. We therefore test the other 
components by imposing a known velocity field, and comparing 
the code's computation of the stress tensor to the analytic solution. 

There are three components remaining to be tested from Equa- 
tionsfToltofBl 

Srr = 2^ (23) 
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(24) 



(25) 
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(b) Angular momentum conservation as a function of r 



Figure 3. Conservation of mass and angular momentum in the FLASH cal- 
culation show in Figure l2l 



of the viscously spreading ring makes all of these terms vanish. We 
therefore adopt a perturbed field: 



v r — vi sin(mi0) + v% sin(fcr) 
i)0 = \/GM,r" 5 + 773 sin(m30 



(26) 
(27) 



The Keplerian velocity field which describes the initial conditions 



These may be easily substituted into Equations 123 1 to l25l to derive 
appropriate analytic expressions. Evaluation of S rr is straightfor- 
ward, since it only involves radial terms. Computing requires 
some care, to ensure that v r is evaluated on a cf> face - a simple 
average is sufficient. Unfortunately, B must be evaluated on two 
different faces (similar to the complication with T r ^, and T^ P ), so 
two separate routines are needed. 

Comparing the code output to the analytic expressions, we 
have found that S rr is accurate to < 1% everywhere except the 
inner boundary (where differences are to be expected). The 
component is also accurate to a similar degree, everywhere except 
where it is close to zero (and relative accuracy becomes meaning- 
less). Performance evaluating B was not quite so good, although 
still satisfactory. Numerical acrobatics similar to those involved in 
Equation ll8l are required to get the derivatives centred correctly. 



6.3 Which components are needed? 

As noted in the introduction, one of the major goals in including 
a physical viscosity into hydrodynamic calculations is to approxi- 
mate the effects of MHD turbulence. Unfortunately, it is not certain 
how best to attain this. If the disc is to accrete, then the T r( f, compo- 
nent of the stress tensor is certainly required. But what of the other 
components? For examp le, in sufficiently v iscous flows, convec- 
tion can be shut down (cf Naravan et all200d) if all components of 
the stress tensor are considered. This is because the convective mo- 
tions are subsonic (in contrast to the supersonic Keplerian shear), 
and hence have very low Reynolds numbers. The effects of includ- 
ing different components of the stress tensor must be considered 
very carefully when performing numerical experiments. 



7 CONCLUSIONS 

In this paper, we have presented an alternative method for including 
viscosity into hydrodynamics codes. It is based on the adjustment 
of velocity fluxes, immediately prior to the transport step. This ap- 
proach ensures excellent conservation of angular momentum, when 
used in polar co-ordinates. We have demonstrated the usefulness 
of this algorithm in the FLASH code. The computational cost is 
low, as compared to the Rieman solver. Additionally, the scheme 
only requires first derivatives, which is numerically desirable. Al- 
though this general approach to solving diffusion equations is not 
new, we are unaware of any contemporary hydrodynamics codes 
which make use of it. In this paper, we have demonstrated the po- 
tential power of the technique. 
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